/*
 heat-spectral lib. solves time-dependent PDEs using SEM-NI method.
 Copyright (C) 2009, Altieri Michele and Stefanoni Andrea
 
 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.
 
 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 GNU General Public License for more details.
 
 You should have received a copy of the GNU General Public License
 along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

#include<ctime>
#include<iostream>
#include<list>

#include "exact.hpp"

using namespace spectral;
using namespace std;

MAKE_L2TIME_ERROR(ExactSolution,Mesh) myL2error;
MAKE_H1TIME_ERROR(ExactSolution,ExactSolutionX,ExactSolutionY,Mesh) myH1error;

int main() {
    // =============================================
    // Problem Parameters
    // =============================================
    // Number of elements along the x-axis
    int nx(4);
    // Number of elements along the y-axis
    int ny(4);
    // Number of nodes
    int nnodes(4);
    // Time step
    double deltat(0.001);
    // Number of time steps
    int ntime(1000);
    // Relaxation parameter
    double theta(0.5);
    // nmax for the newton method
    int newton_nmax((int)1e4);
    // Order of the problem in the week form
    int order(1);
    // Tolerance parameters for the newton method
    double legDelta(1e-15);
    double legEps(1e-15);
    // =============================================

    // =============================================
    // Solving the problem
    // =============================================    
    // Local variables
    double actualtime(0);
    double L2error(0);
    double H1error(0);
    double LinfL2error(0);
    double L2H1error(0);
    double solMin(0);
    double solMax(0);
    char filename[30];
    ofstream l2errorfile;
    ofstream h1errorfile;
    clock_t ck0;
    clock_t ck1;

    // Setting the files for printing the errors
    l2errorfile.open("L2error.dat", std::ios_base::out);
    h1errorfile.open("H1error.dat", std::ios_base::out);
    
    cout << "Initializing SEM problem...\n" << flush;
    ck0=clock();
    Problem<EXACT_PROBLEM> pb(LA,LB,nx,ny,nnodes,deltat,theta,newton_nmax,
			     legDelta,legEps,order);
    pb.initialize();
    pb.printSol("./exactData/exactSol0000.dat");

    ck1=clock();
    cout << "done! ";
    cout << "(In " << double(ck1-ck0)/CLOCKS_PER_SEC << " seconds)\n";
    cout << flush;

    // Computing the range of the initial data
    solMin = pb.range(0);
    solMax = pb.range(1);

    // Evaluating and writing the initial errors
    L2error = sqrt(max(0.,pb.evaluate<myL2error>(0)));
    H1error = sqrt(max(0.,pb.evaluate<myH1error>(0)));

    LinfL2error = L2error;
    L2H1error = deltat/2 * pow(H1error,2);

    l2errorfile << L2error << endl;
    h1errorfile << H1error << endl;
    
    cout << "Solving parabolic problem..." << endl;
    
    ck0=clock();
    for (int t=0;t<ntime;++t) {
	actualtime = deltat*(t+1);

	cout << "Resolution of the time " << actualtime << endl;
	pb.solve();

	if (pb.range(0) < solMin) solMin = pb.range(0);
	if (pb.range(1) > solMax) solMax = pb.range(1);

	sprintf(filename, "%s%04d%s", "./exactData/exactSol", t+1, ".dat");
	pb.printSol(filename);
	
	// Computes the error in L2 norm:
	L2error = sqrt(max(0.,pb.evaluate<myL2error>(actualtime)));
	LinfL2error = max(LinfL2error, L2error);

	l2errorfile << L2error << endl;
	
	// Computes the error in H1 norm:
	H1error = sqrt(max(0.,pb.evaluate<myH1error>(actualtime)));
	L2H1error = L2H1error + (deltat * pow(H1error,2));

	h1errorfile << H1error << endl;
    }
    ck1=clock();
    cout << "done! ";
    cout << "(In " << double(ck1-ck0)/CLOCKS_PER_SEC << " seconds)\n";
    cout << flush;

    // Correction of the L2(H1) norm:
    L2H1error = L2H1error - (deltat/2 * pow(H1error,2));

    l2errorfile.close();
    h1errorfile.close();

    pb.printAviScript("exact", "Heat Problem", ntime, solMin, solMax);
    pb.printGnuplotScript("exact", "Heat Problem", solMin, solMax, true, true);
    pb.printPattern();
    // =============================================

    // =============================================
    // Results output
    // =============================================
    cout << endl;
    cout << "L2(T) Error    : " << L2error << endl;
    cout << "H1(T) Error    : " << H1error << endl;
    cout << "Linf(L2) Error : " << LinfL2error << endl;
    cout << "L2(H1) Error   : " << sqrt(max(0.,L2H1error)) << endl;
    // =============================================
}
